module mod_ode_euler
   use petsc_wrapper
   use mod_ode_solver_abstract
   use stdlib_optval
   use iso_fortran_env, only: dp => real64

   private
   public :: type_ode_euler

   type, extends(type_ode_solver_abstract) :: type_ode_euler

   contains

      procedure :: initialize
      procedure :: destroy
      #: for i in range(1, 5)
      procedure :: step_array${i}$D
      procedure :: substep_array${i}$D
      #: endfor
      procedure :: step_vector
      procedure :: substep_vector

      procedure, private :: time_advance => TimeAdvanceEuler
      #: for i in range(1, 5)
      procedure, private :: time_advance_${i}$d => TimeAdvanceEuler${i}$D
      #: endfor
      procedure, private :: time_advance_vector => TimeAdvanceEulerVector

   end type type_ode_euler

contains

   subroutine initialize(this, scheme)
      class(type_ode_euler), intent(inout) :: this
      character(len=*), intent(in), optional :: scheme
      this%name = "Euler"
      this%scheme = optval(scheme, "explicit")

      select case (this%scheme)
      case ("explicit")
         call petsc_print("Euler scheme: explicit")
         this%order = 1
      case ("implicit")
         call petsc_print("Euler scheme: implicit")
         this%order = 1
      case default
         call petsc_error_print("unknown Euler scheme "//this%scheme//".")
      end select
   end subroutine initialize

   subroutine destroy(this)
      class(type_ode_euler), intent(inout) :: this

      this%A_scalar => null()
      this%A_vector => null()

   end subroutine destroy

   #: for i, dimension in list(zip(range(1, 5), [":", ":, :", ":, :, :", ":, :, :, :"]))
   subroutine step_array${i}$D(this, solution, ctx)
      class(type_ode_euler), intent(inout) :: this
      real(dp), intent(inout), target :: solution(${dimension}$)
      type(c_ptr), intent(inout) :: ctx
      integer :: lrk
      real(dp) :: t

      this%solution_array${i}$_n(${dimension}$) = solution(${dimension}$)
      this%solution_array${i}$ => solution
      call this%solve_rhs(t, solution, ctx)
      call this%time_advance(lrk)

   end subroutine step_array${i}$D

   subroutine substep_array${i}$D(this, lrk, solution, ctx)
      class(type_ode_euler), intent(inout) :: this
      real(dp), intent(inout), target :: solution(${dimension}$)
      type(c_ptr), intent(inout) :: ctx
      integer, intent(in) :: lrk
      real(dp) :: t

      if(lrk == 1) then
         this%solution_array${i}$_n(${dimension}$) = solution(${dimension}$)
      endif
      this%solution_array${i}$ => solution
      call this%solve_rhs(t, solution, ctx)
      call this%time_advance(lrk)

   end subroutine substep_array${i}$D
   #: endfor

   subroutine step_vector(this, solution, ctx)
      class(type_ode_euler), intent(inout) :: this
      class(type_domain_vector), intent(inout), target :: solution
      type(c_ptr), intent(inout) :: ctx
      integer :: lrk
      real(dp) :: t

      if (.not. this%solution_vector_n%allocated) call this%solution_vector_n%duplicate(solution)
      call this%solution_vector_n%copy(solution)
      this%solution_vector => solution
      call this%solve_rhs(t, solution, ctx)
      call this%time_advance(lrk)

   end subroutine step_vector

   subroutine substep_vector(this, lrk, solution, ctx)
      class(type_ode_euler), intent(inout) :: this
      class(type_domain_vector), intent(inout), target :: solution
      type(c_ptr), intent(inout) :: ctx
      integer, intent(in) :: lrk
      real(dp) :: t

      if(lrk == 1) then
         if (.not. this%solution_vector_n%allocated) call this%solution_vector_n%duplicate(solution)
         call this%solution_vector_n%copy(solution)
      endif

      this%solution_vector => solution
      call this%solve_rhs(t, solution, ctx)
      call this%time_advance(lrk)

   end subroutine substep_vector

   subroutine TimeAdvanceEuler(this, lrk)
      class(type_ode_euler), intent(inout) :: this
      integer, intent(in) :: lrk

      select case (this%data_type)
         #: for i in range(1, 5)
      case ("array${i}$D")
         call this%time_advance_${i}$d(lrk)
         #: endfor
      case ("vector")
         call this%time_advance_vector(lrk)
      end select

   end subroutine TimeAdvanceEuler

   #: for i in range(1, 5)
   subroutine TimeAdvanceEuler${i}$D(this, lrk)
      class(type_ode_euler), intent(inout) :: this
      integer, intent(in) :: lrk
      integer :: i, n
      real(dp), allocatable :: A(:, :)

      this%solution_array${i}$ = this%rhs_${i}$*this%dt
      if(this%scheme == "explicit")then
         this%solution_array${i}$ = this%solution_array${i}$_n + this%solution_array${i}$
      elseif(this%scheme == "implicit") then
         if(this%data_type == "array1D") then
            n = size(this%solution_array1)
            if(associated(this%A_scalar))then
               allocate(A, source = this%A_scalar)
               A = this%A_scalar*this%dt 
            else
               allocate(A(n, n)); A = 0.0_dp
            endif
            do i = 1, n
               A(i, i) = A(i, i) - 1.0_dp
            enddo
            call gauss(n, -A, this%solution_array1, this%solution_array1)
            deallocate(A)
         endif
         this%solution_array${i}$ = this%solution_array${i}$_n + this%solution_array${i}$
      endif

   end subroutine TimeAdvanceEuler${i}$D
   #: endfor

! ==========================================================================
!  Ax=b  Solver                                                            =
!  Gauss  Elimination with Maximal Column Pivoting                         =            
!  Gauss  列主元消去法， 求解Ax=b                                          =
!  由于消元过程，计算后 A和b 的值会发生改变                                =
! ========================================================================== 
   subroutine Gauss(n,A,b,x)
      implicit none
      integer:: n,i,j,k,kk
      real(dp):: A(n,n),b(n),x(n)
      real(dp):: fmax,t,f
     do k=1,n-1
  ! ---- find main element ---        
       fmax=abs(a(k,k))
       kk=k
       do i=k+1,n
       if( abs(a(i,k)).gt.fmax )  then
        fmax=abs(a(i,k))
        kk=i
       endif
       enddo
  !c ----- exchange line------         
       do j=k,n
       t=a(k,j)
       a(k,j)=a(kk,j)
       a(kk,j)=t
       enddo
       t=b(k)
       b(k)=b(kk)
       b(kk)=t
  !c    -----Elimination --
      do i=k+1,n
       f=-a(i,k)/a(k,k)
       do j=1,n
       a(i,j)=a(i,j)+f*a(k,j)
       enddo
       b(i)=b(i)+f*b(k)
      enddo
     enddo
  !c  ----- get x ---------       
     x(n)=b(n)/a(n,n)
     do i=n-1,1,-1
     x(i)=b(i)
     do j=i+1,n
     x(i)=x(i)-a(i,j)*x(j)
     enddo
     x(i)=x(i)/a(i,i)
      enddo
  
   end  

   subroutine TimeAdvanceEulerVector(this, lrk)
      class(type_ode_euler), intent(inout) :: this
      integer, intent(in) :: lrk

      call this%solution_vector%copy(this%rhs_vector)
      call this%solution_vector%scale(this%dt)
      if(this%scheme == "explicit") then
         call this%solution_vector%axpy(1.0_dp, this%solution_vector_n)
      elseif(this%scheme == "implicit")then
      endif

      ! call this%solution_vector%axpbypcz(this%alpha(lrk), this%beta(lrk)*this%dt, this%beta(lrk), &
      !                                    this%solution_vector_n, this%rhs_vector)
   end subroutine TimeAdvanceEulerVector

end module mod_ode_euler
